Nuclear size correction to the Lamb shift of one-electron atoms 
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The nuclear size effect on the one-loop self energy and vacuum polarization is evaluated for the 
Is, 2s, 3s, 2pi/ 2 , and 2p 3 / 2 states of hydrogen-like ions. The calculation is performed to all orders 
in the binding nuclear strength parameter Za. Detailed comparison is made with previous all-order 
calculations and calculations based on the expansion in the parameter Za. Extrapolation of the 



all-order numerical results obtained towards Z 
effect on the hydrogen Lamb shift. 
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Introduction 



The distribution of charge of the nucleus influences the 
Dirac energies of atomic systems as well as the quan- 
tum electrodynamical corrections to the energy levels. 
This effect, termed as the nuclear size (NS) correction, 
is important for comparison of theoretical predictions 
with experimental data for the whole range of the nu- 
clear charges Z, from hydrogen (Z = 1) up to uranium 
(Z = 92). The NS corrections to the one- loop self en- 
ergy and vacuum polarization have been previously in- 
vestigated both within the approach based on the expan- 
sion in the nuclear binding strength parameter Za [lHy] 
(where a is the fine structure constant) and within the 
numerical approach that accounts the parameter Za to 
all orders 7, 8]. In the high-Z region, the numerical all- 
order approach provides accurate predictions for the NS 
effect on the radiative corrections. For lower Z, however, 
the NS effect diminishes and becomes increasingly more 
difficult to identify in a numerical calculation. On the 
contrary, the Zct-expansion approach provides accurate 
predictions for the \ow-Z ions and only the qualitative es- 
timates in the high-Z region. A quantitative cross-check 
of the two complementary approaches is not simple and 
has not been accomplished up to now. 

The NS corrections became of particular interest re- 
cently, after the results of the muonic hydrogen Lamb- 
shift experiment were announced @. It turned out that 
the value for the proton charge radius r p deduced from 
the muonic hydrogen differs by five standard deviations 
from the spectroscopic value of r p derived from the hy- 
drogen atom. This unexplained disagreement stimulates 
the scientific community to double-check all contribu- 
tions originating from the nuclear charge distribution, 
both for the muonic and normal atoms. 

The aim of the present investigation is to perform an 
accurate numerical all-order calculation of the NS cor- 
rection to the one-loop self energy and vacuum polariza- 
tion and to make a detailed comparison with the Za- 
expansion results available. 



= 1 provides results for the radiative nuclear size 



I. NS CORRECTION TO DIRAC ENERGY 



The leading-order NS correction to the energy levels of 
a hydrogen-like atom is defined as the difference of the 
corresponding eigenvalues of the Dirac equation with the 
point-Coulomb and the extended-nucleus potentials. The 
two most commonly used models of the nuclear charge 
distribution are the uniformly charged sphere ("sph") 
and the two-parameter Fermi ("Fer") model, 
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where i? S ph = v/5/3i2 is the radius of the sphere with 
the root-mean-square (rms) radius R, c and a are the 
parameters of the Fermi distribution, and po & r e the nor- 
malization factors. The parameter a of the Fermi dis- 
tribution is standardly fixed by a — 2.3/(4 In 3) fm. For 
a given value of the rms radius, the parameter c can be 
determined by the simple approximate formula 
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The relativistic units (to = h - 



1) and the charge 



units a = e /(An) are used in this paper. 



In the calculations performed in this work, we will assume 
that the parameter c of the Fermi distribution is fixed 
exactly by the above formula. 

For the uniformly charged sphere model, the Dirac 
equation can be solved analytically [l(|. In this case, 
the NS correction to the Dirac energy is represented in 
terms of the hypergeometric function. While the exact 
expression is rather cumbersome, simple approximate ex- 
pressions for the NS correction obtained in Ref. [10( are 
highly useful. For the general model of the nuclear charge 
distribution, the NS correction can be easily obtained 
by numerical solution of the Dirac equation. Numerical 
results are conveniently parameterized in terms of the 
function G^(Za, R) 7 whose definition is inspired by the 
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analytic relativistic results (loj . 
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where 7 = -^/l — (Za) 2 and -R sp h is the radius of the 
sphere with the rms radius R. The function Gjq(Za, R) is 
a slowly varying function of Za and R and its numerical 
values are of order of unity. 

The numerical results obtained for the NS correction 
with the Fermi model of the nuclear charge distribution 
are listed in Table [I] The numerical evaluation was per- 
formed by solving the Dirac equation with help of the 
RADIAL package [ll[ and, independently, by using the 
B-spline finite basis set method For calculations in 
the low-Z region, the RADIAL package was upgraded 
into the quadruple arithmetics (about 32 digits). The 
values of the rms charge radii used were taken from the 
compilation 131 for all ions except for uranium; the ura- 
nium rms radius was taken from Ref. [l4j . 



II. NS CORRECTION TO SELF ENERGY 

The one-loop self-energy contribution to the Lamb 
shift is given by a matrix element of the self-energy op- 
erator with the mass renormalization part subtracted, 
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where S(e) = E(e) — 5m and 5m is the mass counterterm. 
The self-energy operator is fl5| 

/oo 
duj D^(uj,x 12 ) 
-OO 

x a fJi G(e — w, xi^X2)a u , 
(7) 

where is the photon propagator, G is the Dirac- 
Coulomb Green function, G(e) = (e — H) , H is the 
Dirac-Coulomb Hamiltonian, and a? = 7°7 AI are the 
Dirac matrices. The nuclear-size self-energy (NSE) cor- 
rection is defined as the difference between the matrix 
elements ^ evaluated with the point-Coulomb potential 
and the potential of the extended-charge nucleus. 

Numerical, all-order (in Za) evaluation of the one-loop 
self-energy correction have been extensively discussed in 
the literature over past decades d, ll6l - [2"oT | , both for the 
case of the point- Coulomb and extended-nucleus poten- 
tials. In the present investigation, we employ the method 
developed in our previous work [2l[ for the case of the 
point nucleus. This method can be immediately extended 
to a general (local and spherically-symmetrical) poten- 
tial, provided that one can calculate the Green function 



of the Dirac equation with this potential. (Beside the 
full Green function, the one-potential Green function is 
also needed in actual calculations.) In the present work, 
we develop an efficient scheme of computation of the 
Dirac Green function for a general potential, which is 
described in Appendix [Al for the full Green function and 
Appendix |B] for the one-potential Green function. 

The main advantage of the method reported in 
Ref. [2l[ is a fast convergence of the partial-wave expan- 
sion of the matrix element ([5]). In the present work, we 
calculate the difference between the point-nucleus and 
extended-nucleus matrix elements. For this difference, 
the partial-wave expansion converges even faster (espe- 
cially, in the low-Z region) than for the self-energy cor- 
rection. Because of this, we were able to significantly 
improve numerical accuracy as compared to results pre- 
viously reported in the literature. 

Numerical results for the NSE correction to the energy 
shift are usually parameterized in the same way as the 
one-loop self-energy itself, 

A£ NSE = - ^nse (Za,R). (8) 

Comparison of the present results with those by Mohr 
and Soff [8j for the homogeneously charged sphere model 
is given in Table |TTJ Numerical results obtained in the 
present work with the Fermi model of the nuclear charge 
distribution are summarized in Table IIIII 

The leading dependence of the NSE correction on R 
and Za can be conveniently factorized out in terms of 
the first-order NS contribution AE-^, 

AE NSE (njl) = AEs(ni/2t) - G NSE (njl) , (9) 

where A£?n is given by Eqs. and (J5J) . An important 
feature of this parametrization [6j, |22j is that it involves 
the full NS correction, rather than only the leading term 
of its Za expansion. With such choice of normalization, 
Gnse is a slowly-varying function of Z and its depen- 
dence on R is more tractable. Note that for the nps/2 
reference state, Eq. (JSJ has AE^(np 1 / 2 ) as a prefactor, 
which was suggested in Ref. [6]. The Za expansion of 
the function Gnse has the form 
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aiog5ji/ 2 In 
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where b = exp[l/(2 7 ) - C - 5/6], 7 = yj\ - (Za) 2 , and 
C is the Euler constant. Known results for the coeffi- 
cients of the expansion are listed in Table IIV1 We note 
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that the logarithmic 022 and 021 terms have not yet ap- 
peared in the literature. It was, however, pointed out 
by Pachucki [23| that such terms are present and that 
the coefficients for the leading logarithms (In 2 (Zee) for 
ns states and hi(Za) for npin states) are the same as 
for the self-energy correction to the hyperfine splitting. 
Values of aoo(?iPi/2) for n > 2 can be found in Ref. 

Comparison of the present numerical data for the func- 
tion Gnse with the Za-expansion results is given in three 
upper graphs of Fig. Q] Note that for ns states, the ratio 
GNSE(Za)/(Za) is plotted. The lower graphs in Fig. Q] 
depict the higher-order remainder (i.e., the contribution 
beyond the known terms of the Za expansion). For ns 
states, the remainder does not approach a finite limit as 
Z — > because it contains ln(Za), as can be seen from 
Eqs. (fTUj) and (fTT|). 

In Fig. [2 the dependence of Gnse on the rms nuclear 
charge radius i? is studied, with the nuclear charge num- 
ber fixed by Z = 92. We find that the i? dependence of 
our numerical data can be well approximated by a three- 
parameter fit that includes In i?, as suggested by the Za 
expansion. More specifically, the following fitting func- 
tions approximate the numerical data for Z = 92 in the 
region R = 3-12 fm with the relative accuracy of better 
than 2 x 10~ 4 (with R expressed in fermi units), 



GnseOR; Is) = -11.8768 + 1.2083 lni? + 0.0191 R, 



(12) 



G N seOR;2s) = -12.2394+ 1.2124 lni? + 0.0220 R, 



(13) 



Gnse 3s) = -11.9131 + 1.2143 lni? + 0.0218 R, 



(14) 

G N SE(i?; 2pi /2 ) = -9.4115 + 1.1759 In i? + 0.0098 i? , 

(15) 

G N SE(i?; 2p 3/2 ) = -1.0145 + 0.0019 In i? + 0.0033 i? . 

(16) 



III. NS CORRECTION TO VACUUM 
POLARIZATION 

The one- loop vacuum-polarization correction to the en- 
ergy levels is usually represented as a sum of two parts, 
the Uehling and the Wichmann-Kroll ones 15]. The 
Uehling part of the vacuum polarization is given by the 
expectation value of the potential 

tt < \ 2a2z r j ' > < >\ 

U\jeh(r) = — / dr r p{r ) 

3mr J 

x [K {2m\r - r'\) - K (2m\r + r'Q] , (17) 



where 

K (x) 



'7^ ^) V^l, MM 



and the nuclear-charge density p is normalized by the 
condition J dr p{r) = 1. The energy shift due to the 



Wichmann-Kroll part of the vacuum polarization can be 
written as @, H3] 



AE- 



WK 



c\ POO POO 

— Re ^ |«| / du / r*dr(gl + f a ) 

x dr'r' (1-7) TrGl 2+ ^y,r'), 

(19) 

where g a and f a are the upper and the lower radial com- 

(24-) 

ponents of the reference-state wave function, G K is the 
radial Dirac Green function that contains two and more 
interactions with the binding field, 

POO 

Gl 2 +) (u,x,y)= / dz z 2 Gf (w, x, z) V(z) 



(20) 



G K is the radial part of the full Dirac Green function, 
Gk "* is the free Dirac Green function, and V(z) is the 
binding potential. We note that Eq. (fl9|) is valid both 
for the point-nucleus and the extended-nucleus binding 
potentials. 

Calculations of the Wichmann-Kroll part of the one- 
loop vacuum polarization were extensively discussed in 
the literature over past decades [7J, |24j-|26( . In the present 
work, we perform calculations of the vacuum polariza- 
tion, evaluating the integrations and the summation over 
k in the order specified by Eqs. (fT9| and (|20|) . Compar- 
ison of the numerical results obtained in this work for 
the Wichmann-Kroll correction with those reported in 
previous calculations [H, [25[ is presented in Table [V] 

The nuclear-size vacuum-polarization (NVP) correc- 
tion is defined as the difference between the one- 
loop vacuum-polarization corrections evaluated with 
the point-Coulomb potential and the potential of the 
extended-charge nucleus. The NVP correction can be 
parameterized in the same way as the one-loop radiative 
corrections, 

A£ NV p - - F NVP (Za, R) . (21) 

The results of our numerical evaluation of the NVP cor- 
rection for the Is, 2s, 3s, 2pi/ 2 , and 2p 3 / 2 states of H- 
like ions are presented in Tabl dVIl The calculation is 
performed for the Fermi model of the nuclear charge dis- 
tribution. It is interesting to note that for the 2p 3 / 2 state 
and high nuclear charges, the correction coming from the 
Wichmann-Kroll part of the vacuum polarization domi- 
nates over the Uehling part. 

The leading dependence of the NVP correction on R 
and Za can be conveniently factorized out in terms of 
the first-order NS contribution A^n BUI, 



AE NVP (njl) = A£ N (ni/ 2 - G NV p(n?7) . (22) 

7T 
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Note that similarly to the NSE correction, for the np^/ 2 
reference state, Eq. (|2"2"| has the first-order NS correction 
for the npi/2 state as a prefactor. 

The Za expansion of the function Gnvp is given by 



Gnvp(«s) = (Za) a w + (Zee)' 



30 



hr 



R 



sph 



+ 02i ln(Za)- 2 + f(Za, R sph ) + 0(1) 



G N vr>(npj) =a 00 + (Za) a w + (Za) 2 

+ 5 j , 1/2 f(Za,R sph )+0(l) 



(23) 

37 -R sP h 



(24) 



where b = exp[l/(2 7 )-C-5/6], 7 = yjl - (Za) 2 , and G 
is the Euler constant. The leading term of the expansion 
for the s states was calculated long ago 0, [1J . All other 
coefficients except 021 were derived in Refs. 0, @. The 
logarithmic 021 term was pointed out by Pachucki [23j j ; 
the value of the coefficient is the same as for the vacuum- 
polarization correction to the hyperfine splitting. The 
results for the expansion coefficients are 

aoa(np 1 / 2 ) = a o(np 3 / 2 ) = - 8 /i5 , (25) 
a 10 (ns) = 3tt/4 , a 10 (np 1/2 ) = 2 ^/72 , a 10 (np 3/2 ) = 5 V 72 > 

(26) 

02i (ns) = 4 /is , (27) 



and 



f(Za, R sph ) = 



1 



3(Za) 2 
2 



-2 ln i? sph - - + 7r tan(7T7) 

3 + 27+ 2^(-l-27) 

7r 3 / 2 (3 + 2 7 )r(7 + f) 

40 sin(27T7)(7 - 1)T(-1 - 2 7 )r( 7 + 3/2) 



x (2i? sph ) 



2(1-7) 



(28) 



The function f(Za, R sp h) has a finite limit as Za — > 0, 
which is 



/(0, -Rsph) — g hi 2 R sph 



4 2 
- - + -C) \nR sph 



+ l -(^- l *C + C 2 -*-\. (29) 
3 V255 5 12 / V ' 

In Fig. [3l we compare the present numerical data for 
the function Gnvp with the Za-expansion results sum- 
marized above. We observe good agreement in all cases; 
the higher-order remainder function exhibits a nearly lin- 
ear dependence on the nuclear charge number. 

In Fig. SI we study the dependence of the function 
Gnvp on the rms nuclear charge radius R, with the nu- 
clear charge number fixed by Z = 92. Similarly to the 



NSE correction, we find that the R dependence of our 
numerical data can be well approximated by a three- 
parameter fit, whose form is suggested by the Za ex- 
pansion. More specifically, the following fitting functions 
approximate the numerical data for Z = 92 in the region 
R = 3-12 fm with the relative accuracy of better than 
2 x 10 -4 (with R expressed in fermi units), 



Gnvp (R; Is) = 15.3607 + 0.3459 In 2 R- 4.4325 InR, 



(30) 



Gnvp (R; 2s) = 15.5292 + 0.3397 In 2 R - 4.4307 ln R , 



(31) 

Gnvp (R] 3s) = 15.4820 + 0.3396 In 2 R - 4.4321 ln R , 

(32) 

G N vpOR;2pi /2 ) = 14.3668 + 0.3673 In 2 R- 4.4346 lni?, 

(33) 

G N vp(-R;2p 3/2 ) = -0.02474- 0.000134 R + 0.000001 R 2 . 

(34) 



IV. RESULTS FOR HYDROGEN 

In this section, we obtain all-order (in Za) results for 
the radiative nuclear size effect to the ground-state Lamb 
shift in hydrogen. This task is complicated by the fact 
that we are not able to perform calculations of the self- 
energy and Wichmann-Kroll parts of the nuclear size ef- 
fect directly for Z = 1. In the absence of a direct cal- 
culation, we perform extrapolation of the numerical data 
obtained for higher values of Z to Z = 1. 

We start with the self energy. The data for the func- 
tion Gnse(^, R) plotted in Fig. [I] is not well suited for 
extrapolation since individual points correspond to dif- 
ferent values of the rms radius. Because of this, we repeat 
our calculations for different nuclear charges and the rms 
radius fixed by R = r p , where r p = 0.8768(69) fm is the 
CODATA value of the proton charge radius [23]. We also 
account for the fact that the Fermi model of the nuclear 
charge distribution is not completely adequate for small 
rms radii; the Gaussian model is employed instead, with 
p(r) — pa exp(— Ar 2 ). The extrapolation is performed 
for the higher-order remainder function, 



£nse = [G N SE(num) - G N SE(ana)]/(Za) , 



(35) 



where GNSE(num) and GNSE(ana) denote the numerical 
and analytical [Eq. (fTTjj) ] values of the Gnse function. In 
our extrapolation, we used 20 points with the nuclear 
charges in the interval Z = 5 — 30 and the same extrapo- 
lation procedure as in Ref. [28j . Our result for hydrogen 
is Gnse( z = 1) = 0.075(25). 

The Uehling part of the NVP correction is calcu- 
lated directly, with the result (for the Gaussian model) 
Gnvp.Uc(Z = 1) = 2.5835 a. The Wichmann-Kroll 
part is a small correction for hydrogen, its leading 
contribution to Gnvp being a constant term of order 
(Za) 2 . Similarly to the NSE correction, we obtain 
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the result for hydrogen by extrapolation. The data to 
be extrapolated is obtained by repeating our calcula- 
tions for Z — 20 — 75, with the rms radius fixed by 
R = r p and the nuclear charge distribution given by 
the Gaussian model. The extrapolation is performed for 
the ratio Gnvp mvk/ {Za} 2 . The result for hydrogen is 
Gnvp,wk(2 = 1) = -9.8(9) a 2 . 

Summarising our calculations of the radiative nuclear 
size effect to the Is Lamb shift in hydrogen, we express 
the results in the same form as in Ref. [23, 



A£ NS e = a 2 Z£ N [-3.1294(80)] , 
AE NV p = a 2 Z£ N [0.8228- 0.0228(23)] 



(36) 
(37) 



where £n = 2/3(Za) 4 R 2 and the first and the sec- 
ond terms in the brackets in Eq. (|3T[) correspond to the 
Uehling and Wichmann-Kroll parts. In order to esti- 
mate the model dependence of our results, we evaluated 
the Uehling part also within the exponential model, with 
p[r) = po exp(— A?'), and found a 0.2% deviation from 
the Gaussian result. 

The numerical constant terms in Eqs. (|36p and (|37[) 
can be compared with the Za-expansion results. For the 
self-energy, the leading-order term of the Za expansion 
is 4 In 2 - 23/ 4 = -2.977, whereas all terms in Eq. (fID|) 
yield the coefficient of —3.153. For the vacuum polar- 
ization, the leading-order term is 3 /4, whereas all terms 
in Eq. fl23j) yield 0.827. We conclude that the higher- 
order corrections increase the leading-order result for the 
radiative nuclear size effect in hydrogen by 4.4%. 



Conclusion 

In the present investigation, we evaluate the nuclear 
size correction to the Lamb shift of the Is, 2s, 3s, 2p 1 / 2 , 
and 2p 3 / 2 states of hydrogen- like atoms. The treatment 
is complete at the one- loop level, i.e., it includes the 
leading-order effect as well as the one-loop radiative cor- 
rections. The total nuclear size correction to the energy 
level is represented, for the ns and npi/ 2 states, as 



A£ns — AE^ 



H — (Gnse + Gnvp) 

TT 



(38) 



and, for the np^/ 2 states, as 
AE- NS (np 3/2 ) =AE- N (np 3/2 ) 

OL 

+ AE N (np 1/2 ) — (G NSE + G NV p) , 

' TT 

(39) 

where A£n is the nuclear size correction to the Dirac 
energy. The all-order numerical values obtained for the 
self-energy and vacuum-polarization functions Gnse and 
Gnvp were compared with results of the Za-expansion 
calculations. Inclusion of the logarithmic term of the rel- 
ative order (Za) 2 hi 2 (Za)~ 2 for ns states was necessary 



in order to achieve agreement between different calcula- 
tions. Extrapolation of the all-order data obtained for 
hydrogen-like ions to Z = 1 provides an all-order result 
for the radiative nuclear size effect on the ground-state 
Lamb shift in hyrogen. The higher-order corrections are 
shown to increase the leading-order result by 4.4%. 
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Appendix A: Dirac Green function for a general 
potential 

In this section we construct the Green function of 
the Dirac equation with the potential V(r) of a general 
form. We assume that the potential V(r) differs from 
the Coulomb one within a finite (inner) region only, i.e., 
there is Tq such that, for r > r , V(r) = —Za/r, with 
Z > 0. For our purposes, it is sufficient to consider the 
potential to be regular at the origin, i.e., rV(r) — >• as 
r — > 0. In the inner region r < ro, the combination 
rV(r) is assumed to be well represented by a piecewise 
cubic polynomial calculated on a sufficiently dense grid. 

The radial Dirac Green function G K (E 1 ri 1 r 2 ) is ex- 
pressed in terms of the two-component solutions of the 
radial Dirac equation regular at the origin (4>^) and the 
infinity (<p^) as follows 

G K (E, n , ra) = - <j%{E, n) <j>f (E, r 2 ) 0(n - ra) 

-4>l{E,r 1 )ct ) f{E 1 r 2 )e{r 2 -T l ). 

(Al) 

The solutions are normalized by the condition that their 
Wronskian is unity (everywhere except for the bound- 
state energies), 



{E, 



-1 

1 



1 



(A2) 



In the present work, we obtain the radial solutions in 
the inner region r < ro by a numerical solution of the 
Dirac equation on a grid, whereas in the outer region 
r > ro, we express them as a combination of the ra- 
dial solutions of the Dirac-Coulomb problem. The reg- 
ular and irregular solutions of the Dirac equation with 
the point-nucleus Coulomb potential will be denoted by 
ip® and ip%*, respectively. They are known analyti cally 
in terms of the Whittaker functions, see e.g., Ref. fig ]. 
(Note the sign difference of the present definition of the 
Green function as compared to that of Ref. [Hj].) In this 
work, the Dirac-Coulomb solutions ^ and tp^ are eval- 
uated by a generalization of the procedure described in 
Ref. [13. 
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TABLE I: Nuclear-size correction to the Dirac energies of H-like ions, in terms of the function Gn(2«, R) defined by Eq. Q. 
Fermi model of the nuclear charge distribution is used. 



z 


n r r 1 

R [fmj 


Is 


2s 


3s 


2P1/2 


5 


2.4059 


1.000 46 


1.000 71 


1.000 23 


1.00173 


8 


2.7013 


1.003 91 


1.004 55 


1.003 33 


1.006 89 


10 


3.0053 


1.006 57 


1.007 58 


1.005 67 


1.011 17 


15 


3.1888 


1.015 66 


1.01793 


1.013 60 


1.025 66 


20 


3.4764 


1.028 67 


1.032 74 


1.024 91 


1.046 42 


26 


3.7371 


1.049 77 


1.056 75 


1.043 18 


1.080 18 


30 


3.9286 


1.06732 


1.076 73 


1.058 28 


1.108 52 


40 


4.2696 


1.124 66 


1.142 02 


1.106 96 


1.202 64 


50 


4.6543 


1.203 59 


1.232 01 


1.17231 


1.337 09 


60 


4.9118 


1.308 62 


1.35181 


1.256 25 


1.524 64 


70 


5.3115 


1.445 02 


1.50715 


1.359 74 


1.784 78 


82 


5.5010 


1.662 15 


1.752 74 


1.51154 


2.236 31 


92 


5.8569 


1.896 75 


2.013 31 


1.655 09 


2.785 73 


100 


5.8570 


2.128 53 


2.263 06 


1.774 54 


3.393 88 



TABLE II: Different calculations of the nuclear-size self-energy correction to the energy levels of H-like ions, in terms of the 
function Fnse(Zci, R) defined by Eq. JSJ, for the homogeneously charged sphere model. 



R [mil 



Is 



2s 



2pi/2 



Ref. 



26 
54 
92 



3.730 
4.826 
5.863 



-0.000 172 122(4) 
-0.000 172(1) 
-0.001 274 870(2) 
-0.001 275(1) 
-0.018 4918(8) 
-0.018 492(1) 



-0.000 173 29(1) 
-0.000 18(1) 
-0.00146169(1) 
-0.001462(1) 
-0.029 089 3(9) 
-0.029 090(2) 



0.000 000 972(2) 
-0.000 00(1) 
-0.000 020 391(2) 
-0.000 021(1) 
-0.002 482 76(8) 
-0.002 483(1) 



[8] 

[8] 



The general calculational scheme is as follows. For a 
given energy argument E, we calculate and store the so- 
lutions ip° and V^f 011 a radial grid {ri}f =1 and then 
obtain the radial Green function for arbitrary radial ar- 
guments by interpolation. Large number of the mesh 
points (N ~ 10 4 ) and a careful choice of the grid allow 
us to minimize losses of accuracy due to interpolation. In 
order to avoid numerical overflow (underflow) when stor- 
ing the regular and irregular solutions for large imaginary 
energies E and large k, all manipulations are performed 
with the "normalized" solutions in which the approxi- 
mate large-r and small-r behaviour is pulled out, 

^ K (E,r)=r-^e- cr cf> K (E,r), (A3) 
^(E,r) = r^e cr ^(E,r), (A4) 

where c = yY— E 2 . Advantages of the normalized solu- 
tions are, first, that they are more suitable for interpo- 
lation than the original solutions and, more importantly, 
that they can be stored and manipulated within the stan- 
dard double precision arithmetics (in the range of k's rel- 
evant for the present investigation, up to \k\ ~ 30). 

In the inner region r < Tq, we calculate the regular 
solution <j) K (or, rather, <f>° K ) by solving the radial Dirac 
equation on a grid as described in the following, starting 
from r = and up to r = r . For r > r , the potential is 
the Coulomb one and the regular solution ifP K is a linear 



combination of the regular and irregular Dirac-Coulomb 
solutions, 

t^{E,r)=a^E,r)+b^{E,r), r>r Q . (A5) 

The coefficients a and b are defined by the condition that 
the two components of 0° are continuous at r — tq. So, 
we determine the coefficients a and b by matching the nu- 
merical and the analytical solutions at r — tq and employ 
the analytical Dirac-Coulomb functions for calculations 
for r > tq. 

The irregular solution </>JJ° in the outer region is just 
the Dirac Coulomb function, 

<P™(E,r) = iP™(E,r), r > r . (A6) 

So, we use the analytical representation for r > tq. For 
smaller r, the irregular solution is calculated by solving 
the Dirac equation on a grid, downward from r = tq to- 
wards r = 0. The normalization of the numerical solution 
is fixed by requiring continuity at the point r = tq. 

We now turn to the problem of solving the Dirac equa- 
tion with the potential V(r) on a grid. In this work, 
we employ the power series solution method, previously 
applied to the Dirac equation by Salvat et al. [ll|. For 
completeness, we give here the description of the method. 
First, let us solve the equation on the interval [r a ,rb] 
with given boundary conditions at r = r a . The situation 
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TABLE III: Nuclear-size self-energy correction to the energy levels of H-like ions, in terms of the function F^se(Zo, R) defined 
by Eq. ((8|. Fermi model of the nuclear charge distribution is used. The nuclear charge rms radii used are listed in Tabled] 



z 


Is 


2s 


3s 


2pi/ 2 


2P3/2 


5 


-0.000 009 962(6) 


-0.000 009 81(2) 








8 


-0.000 020 969(2) 


-0.000 020 53(1) 




0.000 000107(8) 




10 


-0.000 033 242(2) 


-0.000 032 54(1) 


-0.000 032 12(3) 


0.000 000 201(4) 


0.000 000170(2) 


15 


-0.000 060128(2) 


-0.000 059 10(1) 


-0.000 058 21(4) 


0.000 000 412(2) 


0.000 000 345(2) 


20 


-0.000 102 940(2) 


-0.000 102 06(2) 


-0.000 100 38(2) 


0.000 000 699(2) 


0.000 000 587(2) 


26 


-0.000172 537(2) 


-0.000 173 71(3) 


-0.000170 72(3) 


0.000 000 976(2) 


0.000 000 852(2) 


ou 


n nnn o*^s QQftf 

— U.UUU zoo iJoOyZ ) 


— U.UUU z4o Do^O ) 


n nnn o*3Q AKt^V\ 
—u.uuu zoy 4o^o ) 


n nnn nni noi (o\ 

U.UUU UUl \JZ1{Z ) 


n nnn nnn QfiQfo^ 
u.uuu uuu yoy^z ) 


40 


-0.000 481097(2) 


-0.000 51107(2) 


-0.000 502 4(1) 


-0.000 000 726(2) 


0.000 000 396(2) 


50 


-0.000 960 607(2) 


-0.001075 22(2) 


-0.0010578(1) 


-0.000 010 544(2) 


-0.000 003 299(2) 


60 


-0.001830 75(1) 


-0.002 183 28(2) 


-0.002150 3(4) 


-0.000 046 025(4) 


-0.000 014 862(2) 


70 


-0.003 73019(4) 


-0.004 79185(6) 


-0.004 722 0(2) 


-0.000170 450(6) 


-0.000 047078(2) 


82 


-0.008 404 0(2) 


-0.011977 2(2) 


-0.011 796 6(4) 


-0.000 703 91(2) 


-0.000140 462(2) 


92 


-0.018 426 6(7) 


-0.028 986 6(8) 


-0.028 474 0(6) 


-0.002 473 96(6) 


-0.000 337183(2) 


100 


-0.034 060(2) 


-0.058 444(3) 


-0.057172(1) 


-0.006 696 3(4) 


-0.000 602 404(2) 




FIG. 1: (Color online) Nuclear-size self-energy correction, in terms of Gnse defined by Eq. Q, as a function of the nuclear 
charge number Z. The upper graphs depict the ratio Gnse(Z) /{Za) for the ns states and the function Gnse(Z) for the 
npj states, in comparison with the Za-expansion results. The lower graphs show the difference between the all-order and 
Za-expansion results for the function Gnse divided by (Za) 2 . 



r b < r a is allowed and it is assumed that r a > 0. (The radial Dirac equation is (with m = 1) 

special case of r a = will be considered separately.) The k 

G'(r) = — G(r) + (E-V(r) + l)F(r), (A7) 
r 

F'(r) = - F(r) - (E - V(r) - 1) G(r) , (A8) 
r 

where G(r) — rg(r) and F(r) = rf(r) are the upper and 
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FIG. 2: (Color online) Nuclear-size self-energy correction, in terms of GwsE(Za, R) defined by Eq. as a function of the rms 
nuclear charge radius R, for Z = 92. 




FIG. 3: (Color online) Nuclear-size vacuum-polarization correction, in terms of Gnvp{Zo, R) defined by Eq. (|21[) . as a function 
of the nuclear charge number Z. 



lower components of the radial Dirac solution. Introduc- 
ing new variables x = (r — r a )/h and h = ?*6 — ?' a , the 
equation is written as 

(xh + r a )G' x + nhG + UhF-(xh + r a )hF = 0, (A9) 
(xh + r a )F' x - nhF - UhG - (xh + r a )hG = , 

(A10) 

where U — r(V(r) — E). On the given interval, U is 
represented by a cubic polynomial of x, U = X)fe=o u kX k ■ 



The solutions are represented as power series of the form 
G(x) = ]T a n x n , F(x) = ]T b n x n , (All) 

n=0 n=0 

with the coefficients clq and bo determined by the bound- 
ary conditions ao = G(r a ) and bo — F(r a ). The coeffi- 
cients a n and b n are determined by the recurrence rela- 
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FIG. 4: (Color online) Nuclear-size vacuum-polarization correction, in terms of Gnvp{Za, R) denned by Eq. (|21[) , as a function 
of the rms nuclear charge radius R, for Z — 92. 



TABLE IV: Coefficients of the Za expansion of the NSE cor- 
rection in Eq. (fTU|) . 



Term 


State 


Value 


Ref. 


Ooi 


n P:j 


8/9 


[4,5] 




2pi/2 


0.808 879 967(1) 


[4, 5] 




np 3 /2 


ooo(wpi/ 3 ) - 1 


[4,5] 


dlO 


ns 


TV (- 23 /4 + 4 In 2) 


[3, 22, 29] 




npi/2 


7T (379/432 - I6/3 l n 2) 


[6] 




np-i/2 


tv (559/432 - 4 In 2) 




ffllog 


ns, np 1 / 2 


*7s - 15 / 4 


[4, 22] 


0-22 


ns 


" 2 /3 


[23] 


a-21 


nPi/2 


-2(n 2 - l)/n 2 


[23, 30] 


tions (valid for r a ^ 0) 










1 + «)o n -l + ("0 - 


r )6 n -i 




nr a 




•f (ui - h)b n - 2 


+ U2&n-3 + " 3 0n-4 


] , (A12) 



6„ = |Y-n + 1 + n)b n -i + (uq + r Q )a„_i 

nr a 

+ (ui + h)a n -2 + u2«n-3 + usa n -i] . (A13) 

The solutions at the end point are given by the sum of 
the coefficients, 



G{ n ) 



n=0 



(A14) 



In the numerical evaluation, the recurrence relations are 
applied upwards until either the desired precision or the 
upper limit for n (typically, n ma x = 30) is reached. In 
the latter case, the interval is subdivided into two parts 
and the procedure is repeated until the desired accuracy 
is attained. This simple approach allows one to solve the 
equation with accuracy close to the machine precision. 

Now, we consider the special case of r a = 0. In this 
case, the solutions are represented by the power expan- 



sion of the form 
G(x) = x s a ^ n > F ( x ) = xS+t £ b ^ n > ( Al5 ) 



71=0 



n=0 



where the parameters s and t are determined from the 
Dirac equation. For k < 0, we have (for the regular 
potentials considered here) s = \k\ and t = 1. The series 
start with 



a = 1 , b 



l + 2\K\ 

The recursion relations take the form 



(A16) 



na n = (ui - h)b n -2 + u 2 b n - 3 + u 3 o„_4 , (A17) 
(2|k| + n+ l)b n = (h + U\)a n + u 2 a n -i + u 3 a„_ 2 • 

(A18) 

For k > 0, one gets s = n + 1 and t = — 1. The series 
start with 



h-ui 

an = , On = 1 . 

1 + 2k ' 
whereas the recursion relations are 



(A19) 



(2k + n + l)a„ = (h - ui)6„ - u 2 b n -i - u 3 b n - 2 , 

(A20) 

nb n = (ui + h)a n - 2 + it 2 a„_ 3 + u 3 a„_ 4 . (A21) 

Appendix B: One-potential Dirac Green function for 
a general potential 

For the evaluation of the self-energy correction, the 
one-potential Dirac Green function G^ is needed. Its 
radial part is defined as 

/>oo 

G ( £ ) {E,r l ,r 2 )= / dzz 2 V(z)G^(E,r 1 ,z)G^(E,z,r 2 ). 
Jo 

(Bl) 
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TABLE V: Different calculations of the Wichmann-Kroll vacuum-polarization correction [Eq. (|19[) ] to the energy levels of H-like 
ions, in units AE/[(a/w) (Za) 4 /n 3 ]. 



z 


R [fm] 


Model 


Is 


2s 


2Pl/2 


2^3/2 


Ref. 


36 


4.230 


sphere 


0.002 74199(6) 


0.002 833 98(3) 


0.000 083 36(4) 


0.000 027 95(4) 








sphere 


0.002 741 8 


0.002 833 7 


0.000 083 4 


0.000 028 


[25] 






shell 


0.002 7 


0.002 8 


0.000 1 


0.000 


[15] 


54 


4.826 


sphere 


0.005 9212(2) 


0.006 423 4(2) 


0.000 454 79(2) 


0.00011418(5) 






sphere 


0.005 921 1 


0.006 4231 


0.000 454 8 


0.000114 2 


[25] 






shell 


0.005 9 


0.006 4 


0.000 4 


0.000 1 


[15] 


92 


5.860 


sphere 


0.020 679 2(5) 


0.027 252 0(8) 


0.006 8241(4) 


0.000 749 49(6) 






sphere 


0.020 678 9 


0.027 2515 


0.006 824 


0.000 749 4 


[25] 






shell 


0.020 6 


0.027 2 


0.006 8 


0.000 7 


[15] 



TABLE VI: Nuclear-size vacuum-polarization correction to the energy levels of H-like ions, in terms of the function F^yp(Za, R) 
defined by Eq. (|2ip . Fermi model of the nuclear charge distribution is used. The nuclear charge rms radii used are listed in 
Table[I] For a given Z, the upper line ( "Ue" ) corresponds to the Uehling part and the lower line ( "WK" ) , to the Wichmann-Kroll 
part. 



Z 




Is 


2s 


3s 


2P1/2 


2^3/2 


15 


Ue 


0.000 024 856 


0.000 024 968 


0.000 024 921 


0.000 000 020 


-0.000 000 016 


20 


Ue 


0.000 047 62 


0.000 048 21 


0.000 048 12 


0.000 000102 


-0.000 000 034 




WK 


-0.000 006 4(2) 


-0.000 006(1) 


-0.000 006(3) 






26 


Ue 


0.000 089 44 


0.000 09172 


0.000 091 60 


0.000 000402 


-0.000 000 064 




WK 


-0.000 012 58(1) 


-0.000 012 9(1) 


-0.000 013(1) 






30 


Ue 


0.000 131 907(2) 


0.000 136 725(2) 


0.000 136 601(2) 


0.000 000 865 


-0.000 000 092 




WK 


-0.000 018 83(3) 


-0.000 019 55(9) 


-0.000 019 6(7) 






40 


Ue 


0.000 304 304(4) 


0.000 326 352(4) 


0.000 326 510(4) 


0.000 004 205 


-0.000 000 188 




WK 


-0.000 043 43(4) 


-0.000 046 5(1) 


-0.000 046 6(9) 


-0.000 000 9(1) 




50 


Ue 


0.000 674 503(2) 


0.000 756 416(2) 


0.000 758 099(2) 


0.000 016 672 


-0.000 000 342 




WK 


-0.000 092 38(7) 


-0.000 102 71(7) 


-0.000 103 07(9) 


-0.000 003 37(2) 


-0.000 000 20(2) 


60 


Ue 


0.001410 95(1) 


0.001673 01(1) 


0.001680 03(1) 


0.000 05740 


-0.000 000 546 




WK 


-0.000 1816(2) 


-0.000 2118(3) 


-0.000 212 8(3) 


-0.000 01046(2) 


-0.000 000 49(2) 


70 


Ue 


0.003 100 32(1) 


0.003 932 42(2) 


0.003 956 23(2) 


0.000 197 79 


-0.000 000 872 




WK 


-0.000 367 7(4) 


-0.000 454 6(5) 


-0.000 456 8(5) 


-0.000 032 06(4) 


-0.000 001201(4) 


82 


Ue 


0.007 710 07(4) 


0.010 779 92(6) 


0.010 863 65(6) 


0.000 822 122(5) 


-0.000 001310 




WK 


-0.000 8215(7) 


-0.001 106 2(9) 


-0.0011114(9) 


-0.000114 34(9) 


-0.000 003 021(4) 


92 


Ue 


0.018 230 65 


0.028 056 439(2) 


0.028 275 306(4) 


0.002 970 972 


-0.000 001923 




WK 


-0.001762 6(8) 


-0.002 587(1) 


-0.002 594(1) 


-0.000 3610(2) 


-0.000 006 750(3) 


100 


Ue 


0.036 429 910(6) 


0.061 165 53(1) 


0.06154911(1) 


0.008 427011(2) 


-0.000 002 344 




WK 


-0.003 259(2) 


-0.005 184(4) 


-0.005 183(4) 


-0.000 915 4(6) 


-0.000 012100(6) 



where G^ is the free Dirac Green function. Substituting 
the representation (jAll) for G«P into (|B1I) and introduc- 
ing the integral functions 

J°°(r) = f dzz 2 V{z) V f{z) V l{z), (B2) 
Jo 

J?(r)= f dzz 2 V{z) V f{z)^ K {z), (B3) 
Jo 

J:(r)= dzz 2 V(z)^ T (z)^(z), (B4) 

J r 

where ^k ^ Upk°^ denote the regular (irregular) so- 
lutions of the free Dirac equation, we write the one- 



potential Dirac Green function for r\ < r 2 as 

G£\E,r 1 ,r 2 ) = $°(r 1 )vf(r 2 ) + <p° K (r 1 )$f(r 2 ), 

(B5) 

where 

d>°(r) = p~(r) J°°(r) - y£(r) Jf (r) , (B6) 
<TO = <P? (r) Jl° (r) + <p° K (r) (r) . (B7) 

For r*i > r 2 , the one-potential Green function is obtained 
by the symmetry condition, 

Gi 1 HE,r 1 ,r 2 ) = G^ T (E,r 2 ,r 1 ). (B8) 
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Analogously to the approach used for the full Dirac 
Green function, we store the functions ip K and <£> K on a 
radial grid {rj} i=1 and obtain the one-potential Green 
function by interpolation. The integral functions J K are 
evaluated by numerical integration with help of Gauss- 
Legendre quadratures. The integration interval (0, oo) 
is breaked up at the position of the mesh points rj, so 
that only one integral over (0, oo) needs to be evaluated 



for a given value of E. Analogously to the case of the 
full Green function, all manipulations with the regular 
and irregular solutions are carried out after normalizing 
them according to Eqs. (|A3I) and (|A4I) . in order to pre- 
vent numerical overflow and underflow. Similar method 
of computation of the one-potential Green function was 
used long ago by M. Gyulassy in his evaluation of the 
vacuum-polarization [31j . 
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